EDA¶
In [1]:
import pandas as pd
import numpy as np
In [2]:
import os
os.getcwd()
Out[2]:
In [3]:
os.listdir('data')
Out[3]:
In [4]:
all_data = pd.read_csv('data/train.csv')
- Cat_features: label encode
- city: maybe cardinality handler
In [5]:
all_data.info()
All features profiling¶
In [6]:
# Import pandas_profiling - an useful library for EDA
from pandas_profiling import ProfileReport
# Generate profiling for all features
profile = ProfileReport(all_data, title='Profiling Report', explorative = True)
profile
Out[6]:
Compare feature distribution of the two classes¶
In [7]:
# Since the field `enrollee_id` is 100% distinct and have no missing, it can be used as index
all_data.set_index('Id', inplace = True)
In [8]:
cat_features = []
num_features = []
data_types = all_data.drop('Risk_Flag', axis = 1).dtypes
for feature in data_types.index:
if data_types[feature] == 'object':
cat_features.append(feature)
else:
num_features.append(feature)
In [9]:
cat_features
Out[9]:
In [10]:
num_features
Out[10]:
In [11]:
import matplotlib.pyplot as plt
Numerical features distribution¶
- Each numerical features of each class will be visualized in a histogram
- Since the number of classes are imbalanced, the distributions are percentage instead of absolute count
In [12]:
eda_data = {
'defaulted': all_data[all_data['Risk_Flag'] == 1],
'not_defaulted': all_data[all_data['Risk_Flag'] == 0]
}
weights = dict()
for class_ in eda_data:
weights[class_] = len(eda_data[class_].index)
weights[class_] = np.ones(weights[class_]) / weights[class_]
In [13]:
for col in num_features:
fig, ax = plt.subplots(1, 2)
idx = 0
for class_ in eda_data:
eda_data[class_][col].hist(ax = ax[idx], weights = weights[class_])
ax[idx].set_title(class_)
idx += 1
fig.suptitle(col)
plt.tight_layout()
plt.show()
Striking point(s):
- The distributions of the numerical features of both classes very similar
- The
not_defaultedusers'Experiencehas a higher percentage of worker with 20 years experience
Categorical features (top up to 10 values)¶
In [14]:
for col in cat_features:
fig, ax = plt.subplots(1, 2)
idx = 0
for class_ in eda_data:
eda_data[class_][col].value_counts(normalize = True)[:10].plot.bar(ax = ax[idx])
ax[idx].set_title(class_)
idx += 1
fig.suptitle(col)
plt.tight_layout()
plt.show()
Striking point(s):
- The top
Profession,STATEandCITYof the two classes are very different
Model¶
Estimator¶
In [15]:
from catboost import CatBoostClassifier
In [16]:
class CustomizedCatBoostClassifier(CatBoostClassifier):
def fit(self, X, y, *args, **kwargs):
X_train, X_eval, y_train, y_eval = train_test_split(X,
y,
test_size = 0.3)
super().fit(
X_train,
y_train,
eval_set = (X_eval, y_eval),
*args,
**kwargs
)
return self
Preprocessor¶
Since some features has very high cardinality (e.g. CITY has more than 300 values), a CardinalityHandler which restricts the cardinality is useful to speed the training up
In [17]:
from sklearn import set_config
set_config(display='diagram')
In [18]:
from sklearn.base import TransformerMixin, BaseEstimator
class CardinalityHandler(BaseEstimator, TransformerMixin):
def __init__(self, columns, cardinality):
if len(columns) != len(cardinality):
raise Exception("The len of the list of columns does not match the len of the cardinality")
self.columns = columns
self.cardinality = cardinality
self.top_values = dict()
def fit(self, X, y = None):
for i in range(len(self.columns)):
self.top_values[self.columns[i]] = list(X[self.columns[i]].value_counts().index)[:self.cardinality[i]]
if 'unk' not in self.top_values[self.columns[i]]:
self.top_values[self.columns[i]].append('unk')
if 'capped_value' not in self.top_values[self.columns[i]]:
self.top_values[self.columns[i]].append('capped_value')
return self
def transform(self, X, y = None, copy = False):
if copy:
tranform_data = X.copy()
else:
transform_data = X
for col, top_val in self.top_values.items():
transform_data.loc[:, col] = transform_data[col].apply(lambda x: x if x in top_val else 'capped_value')
return transform_data
In [19]:
from sklearn.impute import SimpleImputer
class CustomizedSimpleImputer(SimpleImputer):
def transform(self, X):
data = super().transform(X)
return pd.DataFrame(data, columns = X.columns, index = X.index)
In [20]:
from sklearn.compose import ColumnTransformer
class CustomizedColumnTransformer(ColumnTransformer):
def _hstack(self, Xs):
return pd.concat(Xs, axis = 1)
In [21]:
from sklearn.pipeline import FeatureUnion, Pipeline
class CustomizedFeatureUninon(FeatureUnion):
def _hstack(self, Xs):
cols = [X.columns.tolist() for X in Xs]
cols = np.hstack(cols)
return pd.DataFrame(super()._hstack(Xs), columns = cols)
In [22]:
cat_imputer = CustomizedSimpleImputer(
missing_values = np.nan,
strategy = 'constant',
fill_value = 'unk', # short for unknown
)
columns = ['CITY']
cardinality = [100]
cardinality_handler = CardinalityHandler(columns, cardinality)
cat_handler = Pipeline(steps = [
('imputer', cat_imputer),
('cardinalityHandler', cardinality_handler)
],
verbose = False)
num_imputer = CustomizedSimpleImputer(
missing_values = np.nan,
strategy = 'most_frequent',
)
preprocessor = CustomizedColumnTransformer([
('catHandler', cat_handler, cat_features),
('numHandler', num_imputer, num_features)
],
n_jobs = -1)
In [23]:
fu = CustomizedFeatureUninon([
('featureUnion', preprocessor)
])
In [24]:
pl = Pipeline(steps = [
('preprocessor', fu),
('estimator', CustomizedCatBoostClassifier(cat_features = cat_features,
use_best_model = True,
verbose = True,
metric_period = 100,
early_stopping_rounds = 25,
auto_class_weights = 'Balanced'
))
],
verbose = True)
In [25]:
pl
Out[25]:
Training¶
Split data¶
In [26]:
from sklearn.model_selection import train_test_split
In [27]:
train_data, test_data = train_test_split(all_data, test_size = 0.2)
In [28]:
train_data.head()
Out[28]:
In [29]:
label = train_data['Risk_Flag'].copy()
train_data = train_data.drop('Risk_Flag', axis = 1)
test_label = test_data['Risk_Flag'].copy()
test_data = test_data.drop('Risk_Flag', axis = 1)
In [30]:
train_data.shape
Out[30]:
In [31]:
label.value_counts(normalize = True)
Out[31]:
In [32]:
test_data.shape
Out[32]:
In [33]:
test_label.value_counts(normalize = True)
Out[33]:
Oversampling¶
The data is extremely unbalanced, therefore, oversampling is necessary
In [34]:
cat_col_loc = []
for col in cat_features:
cat_col_loc.append(train_data.columns.get_loc(col))
In [35]:
from imblearn.over_sampling import SMOTENC
X_resampled, y_resampled = SMOTENC(categorical_features = cat_col_loc).fit_resample(train_data, label)
Fit¶
In [36]:
pl.fit(X_resampled, y_resampled)
Out[36]:
Evaluate¶
In [37]:
y_pred = pl.predict(test_data)
In [38]:
from sklearn.metrics import classification_report, roc_auc_score, plot_confusion_matrix
Classification report¶
In [39]:
print(classification_report(test_label, y_pred))
AUC¶
Raw ROC_AUC score¶
In [40]:
roc_auc_score(test_label, y_pred)
Out[40]:
Classes AUC Curve and score¶
In [41]:
from sklearn.metrics import roc_auc_score, roc_curve, auc
from itertools import cycle
def auc_report(model, X, y_true, fig, ax):
classes = model.classes_
y_pred_classes = model.predict_proba(X)
y_pred = model.predict(X)
n_classes = len(classes)
lw = 2
for i in range(len(classes)):
print(f"""{classes[i]}: {roc_auc_score(y_true=(y_true==classes[i]).astype(int), y_score=y_pred_classes[:,i])}""")
# Compute ROC curve and ROC area for each class
fpr = dict()
tpr = dict()
roc_auc = dict()
for i in range(n_classes):
fpr[i], tpr[i], _ = roc_curve(y_true=(y_true==classes[i]).astype(int), y_score=y_pred_classes[:,i])
roc_auc[i] = auc(fpr[i], tpr[i])
all_fpr = np.unique(np.concatenate([fpr[i] for i in range(len(classes))]))
# Then interpolate all ROC curves at this points
mean_tpr = np.zeros_like(all_fpr)
for i in range(len(classes)):
mean_tpr += np.interp(all_fpr, fpr[i], tpr[i])
# Finally average it and compute AUC
mean_tpr /= n_classes
fpr["macro"] = all_fpr
tpr["macro"] = mean_tpr
roc_auc["macro"] = auc(fpr["macro"], tpr["macro"])
# Plot all ROC curves
ax.plot(fpr["macro"], tpr["macro"],
label='macro-average ROC curve (area = {0:0.2f})'
''.format(roc_auc["macro"]),
color='navy', linestyle=':', linewidth=4)
colors = cycle(['aqua', 'darkorange', 'cornflowerblue'])
for i, color in zip(range(n_classes), colors):
plt.plot(fpr[i], tpr[i], color=color, lw=lw,
label='ROC curve of class {0} (area = {1:0.2f})'
''.format(classes[i], roc_auc[i]))
ax.plot([0, 1], [0, 1], 'k--', lw=lw)
ax.set_xlim([0.0, 1.0])
ax.set_ylim([0.0, 1.05])
ax.set_xlabel('False Positive Rate')
ax.set_ylabel('True Positive Rate')
fig.suptitle('Some extension of Receiver operating characteristic to multi-class')
ax.legend(loc="lower right")
# plt.show()
output = {
'macro_auc': roc_auc["macro"],
# f'auc_{classes[0]}': roc_auc[0],
# f'auc_{classes[1]}': roc_auc[1]
}
for i in range(len(classes)):
output[f'auc_{classes[i]}'] = roc_auc[i]
return output
In [42]:
fig, ax = plt.subplots()
# y_pred_classes = pl.predict_proba(test_data.drop('Risk_Flag', axis = 1))
auc_report(pl, test_data, test_label, fig, ax)
Out[42]:
Confusion matrix¶
In [43]:
IC = type('IdentityClassifier', (), {"predict": lambda i : i, "_estimator_type": "classifier"})
plot_confusion_matrix(IC,
y_pred,
test_label,
normalize='true',
values_format='.2%',
cmap = 'Blues'
);
Model explainability¶
Feature importance¶
In [44]:
fi_df = pd.DataFrame({
'feature': train_data.columns,
'feature_importance':pl.steps[-1][1].get_feature_importance()
})
In [45]:
fi_df.set_index('feature').sort_values(by = 'feature_importance').plot.barh()
Out[45]:
Shap value¶
In [46]:
import shap
In [47]:
X = pl.steps[0][1].fit_transform(test_data)
shap.initjs()
explainer = shap.TreeExplainer(pl.steps[-1][1])
shap_values = explainer.shap_values(X)
shap.summary_plot(shap_values,
X.to_numpy(),
feature_names = X.columns
)
Note: The points of categorical features in the graph does not have any color because those features does not have a hierachy (low, high)
Not Defaulted users¶
In [48]:
X = pl.steps[0][1].fit_transform(test_data.loc[test_label[test_label == 0].index])
shap.initjs()
explainer = shap.TreeExplainer(pl.steps[-1][1])
shap_values = explainer.shap_values(X)
shap.summary_plot(shap_values,
X.to_numpy(),
feature_names = X.columns,
# plot_type = 'bar'
)
Defaulted users¶
In [49]:
X = pl.steps[0][1].fit_transform(test_data.loc[test_label[test_label == 1].index])
shap.initjs()
explainer = shap.TreeExplainer(pl.steps[-1][1])
shap_values = explainer.shap_values(X)
shap.summary_plot(shap_values,
X.to_numpy(),
feature_names = X.columns,
# plot_type = 'bar'
)
Error Analysis¶
- Investigate the pattern of the wrongly classified users by looking at their top features
In [50]:
ea_df = pd.DataFrame({
'y_pred': y_pred,
'y_true': test_label
})
In [51]:
ea_df['true'] = ea_df['y_pred'] == ea_df['y_true']
In [52]:
ea_df = pd.concat([ea_df, test_data], axis = 1)
Wrongly classified as Defaulted¶
In [53]:
ea_df[(ea_df['true'] == False) & (ea_df['y_true'] == 1)]['Profession'].value_counts()[:10].plot.barh()
Out[53]:
In [54]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 1)]['Married/Single'].value_counts()[:10].plot.bar()
Out[54]:
In [55]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 1)]['CURRENT_JOB_YRS'].hist()
Out[55]:
In [56]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 1)]['Experience'].hist()
Out[56]:
Wrongly classified as Not Defaulted¶
In [57]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 0)]['Profession'].value_counts()[:10].plot.barh()
Out[57]:
In [58]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 0)]['Married/Single'].value_counts()[:10].plot.bar()
Out[58]:
In [59]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 0)]['CURRENT_JOB_YRS'].hist()
Out[59]:
In [60]:
ea_df[(ea_df['true'] == False) & (ea_df['y_pred'] == 0)]['Experience'].hist()
Out[60]:
Striking point(s):
- Most users'
Professionare not the popular one of neither classes
Conclusion¶
- Overall, the model performance is acceptable (ROC_AUC score, accuracy, each class' auc)
- The precision of
defaultedclass is low, but it is acceptable because the number of this class is very small compared to the others. Therefore, a small percentage of wrongly classified users in thenot defaultedclass can lead to a low precision indefaulted - The model explainability and feature importance are reasonable
What to do next¶
- Do parameter tuning to improve the performance
- Try different preprocessor (e.g. encode the categorical values to decrease the size of the model)
- Try different tree-based models (LightGBM, XGBoost)
Predict¶
In [61]:
predict_data = pd.read_csv('data/predict_data.csv')
In [62]:
predict_data.head()
Out[62]:
In [63]:
predict_data.set_index('ID', inplace = True)
prediction = pl.predict(predict_data)
In [64]:
prediction
Out[64]:
In [65]:
predict_data['prediction'] = prediction
predict_data[['prediction']].to_csv('data/prediction.csv')